* Code to calculate local projection impulse responses with US IRI data
local folder "myfolder"   //or "D:"

if "`folder'"=="myfolder" {
cd "D:\research\D Karadi-SchoenleWursten-Selection\code\Submission\US_CPI\code"
}
* log using "../data/US_LP.log", replace

global path_graphs = "../text/figures/"

use "../data/aggregateData.dta", clear
rename month dm
tsset dm, monthly
foreach var of varlist ext_CPI_FAH_SA_inf {
	gen d`var'=100*D.`var'
}
* drop categoryID
sort dm
save "../data/IRi_data.dta", replace

* Load us data, prepare it for analysis
import delimited using "../data/usdata.csv", delimiters(",") varnames(1) clear
gen dm=ym(year,month)
format dm %tm
order dm, first
tsset dm, monthly
keep dm ff4_hf sp500_hf gss1_hf signnegm_ff4sp500 signposm_ff4sp500 ///
		 pmnegm_ff4sp500 pmposm_ff4sp500 pmnegm_gss1sp500 pmposm_gss1sp500 ///
		 ebpnew ip sp500 baa_spread_moodys_m mortg_spread_m ff gs1 gs2 cpi cpi_core unrate ///
		 ppi ppi_fg ppi_ind ppi_fg_core
		 
* Create log-differences and annualize aggregate variables
foreach var of varlist ip cpi cpi_core unrate ppi* {
	gen l`var' = 100*log(`var')
	gen dl`var' = D.l`var'
}

* Annualize ip and hicp
foreach var of varlist ff gs1 gs2 ebpnew baa_spread_moodys_m mortg_spread_m sp500 {
	gen d`var'=D.`var'
}

** Seasonal adjustment
gen month_nr = month(dofm(dm))
tab month_nr, gen(month_)
local sample "year(dofm(dm))>=1979 & year(dofm(dm))<=2018"

foreach var of varlist dlppi {
	gen `var'_sa= .
	qui regress `var' month_1-month_12 if `sample', noconst
	qui replace `var'_sa=(_b[month_1]+_b[month_2]+_b[month_3]+_b[month_4]+_b[month_5]+_b[month_6]+_b[month_7]+_b[month_8]+_b[month_9]+_b[month_10]+_b[month_11]+_b[month_12])/12 if `sample'		// Add average coefficient
	qui replace `var'_sa=`var'+`var'_sa -_b[month_1]*month_1-_b[month_2]*month_2-_b[month_3]*month_3-_b[month_4]*month_4-_b[month_5]*month_5-_b[month_6]*month_6-_b[month_7]*month_7-_b[month_8]*month_8-_b[month_9]*month_9-_b[month_10]*month_10-_b[month_11]*month_11-_b[month_12]*month_12 if `sample'	// Subtract seasonal component
	gen d`var'_sa = D.`var'_sa
}
	
drop month_*

gen lppi_sa = sum(dlppi_sa)
replace lppi_sa =. if lppi_sa==0

merge 1:1 dm using "../data/IRi_data.dta"
drop _merge

tsset dm, monthly

//Look at some moments
local varlist "dlcpi_core dlcpi dlppi_sa dext_CPI_FAH_SA_inf"
foreach var of varlist `varlist' {
	gen `var'_ann=12*`var'
}
su dlcpi_core_ann dlcpi_ann dlppi_sa_ann dext_CPI_FAH_SA_inf_ann if year(dofm(dm))>=2001 & year(dofm(dm))<=2012


* Set some parameters
local NImpLen = "24"

foreach pre in DL DLu DLd DP DPu DPd {
	gen m`pre'_price_ref_h0_r_SA = `pre'_price_ref_v1_r_SA
}

foreach h of numlist 0/`NImpLen' {
		gen mSize_price_ref_h`h'_r_SA = mDL_price_ref_h`h'_r_SA/mDP_price_ref_h`h'_r_SA
		gen mSizeU_price_ref_h`h'_r_SA = mDLu_price_ref_h`h'_r_SA/mDPu_price_ref_h`h'_r_SA
		gen mSizeD_price_ref_h`h'_r_SA = mDLd_price_ref_h`h'_r_SA/mDPd_price_ref_h`h'_r_SA
}

//turn growth measures into percent terms
foreach var of varlist *DL* *DP* *Size* {
	replace `var'=100*`var'
	gen d`var'=D.`var'
}

* local price_list = "DL_price_raw_v1_r_SA DL_price_ref_v1_r_SA DL_price_reg5_v1_r_SA DL_RTR_SoD_v1_r_new_SA DL_Sref_v1_r_SA dgs1 dff"
* local price_label_list = "Posted_prices Reference_prices Regular_prices Real_expenditures Sales_prices_(ref) 1_year_rate Federal_funds_rate"
* foreach price_var in "`price_list'" {
* 	local price = "`price_var'"
*}

// price 	instrument controls											level_or_cumulative		sample
//"debpnew" "ebpnew" "lip lcpi_core gs1 l(1/12).(lip lcpi_core gs1 ebpnew)" "level" "year(dofm(dm))>=2001 & year(dofm(dm))<=2012"  //Figure 6a
//"dgs1" "ebpnew" "lip lcpi_core gs1 l(1/12).(lip lcpi_core gs1 ebpnew)" "level" "year(dofm(dm))>=2001 & year(dofm(dm))<=2012"  //Figure 6b
//"dlip" "ebpnew" "lip lcpi_core gs1 l(1/12).(lip lcpi_core gs1 ebpnew)" "level" "year(dofm(dm))>=2001 & year(dofm(dm))<=2012"  //Figure 6c
//"dlcpi_core" "ebpnew" "lip lcpi_core gs1 l(1/12).(lip lcpi_core gs1 ebpnew)" "level" "year(dofm(dm))>=2001 & year(dofm(dm))<=2012" //Figure 6d
//"dext_CPI_FAH_SA_inf" "ebpnew" "lip lcpi_core gs1 l(1/12).(lip lcpi_core gs1 ebpnew)" "level" "year(dofm(dm))>=2001 & year(dofm(dm))<=2012" //Figure 6e
//"DL_price_ref_v1_r_SA" "ebpnew" "lip lcpi_core gs1 l(1/12).(lip lcpi_core gs1 ebpnew)" "level" "year(dofm(dm))>=2001 & year(dofm(dm))<=2012" //Figure 6e

//"debpnew" "ebpnew" "lip lcpi_core gs1 "lip lppi_sa gs1 l(1/12).(lip lppi_sa gs1 ebpnew)" "level" "year(dofm(dm))>=1985 & dm<=tm(2015,12)" //Figure 9a
//"dgs1" "ebpnew" "lip lcpi_core gs1 "lip lppi_sa gs1 l(1/12).(lip lppi_sa gs1 ebpnew)" "level" "year(dofm(dm))>=1985 & dm<=tm(2015,12)" //Figure 9b
//"dlppi_sa" "ebpnew" "lip lcpi_core gs1 "lip lppi_sa gs1 l(1/12).(lip lppi_sa gs1 ebpnew)" "level" "year(dofm(dm))>=1985 & dm<=tm(2015,12)" //Figure 9c
//"dlip" "ebpnew" "lip lcpi_core gs1 "lip lppi_sa gs1 l(1/12).(lip lppi_sa gs1 ebpnew)" "level" "year(dofm(dm))>=1985 & dm<=tm(2015,12)" //Figure 9d


local level_or_cumulative = "level"  //"level": needs to cumulate; "cumulative": already cumulated
local price =  "debpnew"  //"dlip"  //"dlppi_sa"  //"DL_price_ref_v1_r_SA"  //"dlcpi_core"  //"dlip"  
local price_label = "Excess bond premium" //"Industrial production"  //"Producer price index"   //"Reference price inflation" //"Food at home"  //"CPI (core)" // "1-year rate"
local instrument = "ebpnew"     //"pmnegm_ff4sp500"  
local instrument_label "ebp" //"Poor man's monetary policy shock" 
local controls = "lip lcpi_core gs1 l(1/12).(lip lcpi_core gs1 ebpnew)"    // "lip lppi_sa gs1 l(1/12).(lip lppi_sa gs1 ebpnew)"
local controls_label "CPI (log), IP (log), 1-year Treasury, ebp (all with 12 lags)"  //"PPI (log), IP (log), 1-year Treasury, ebp (all with 12 lags)"  
local controls_name "cpi_coreipgs1ebp112" //"ppiipgs1ebp112" 

global balanceSample = "1"
local sample = "year(dofm(dm))>=2001 & year(dofm(dm))<=2012"  // "year(dofm(dm))>=1985 & dm<=tm(2015,12)"  //
local sample_str = "1985_2015"  //"2001_2012" //

* Prep results matrix
mat results = J(`NImpLen'+1, 3, .)
mat colnames results = horizon b se

* Balance sample
if "$balanceSample" == "1" {
	if "`level_or_cumulative'" == "level" {
		reg `price' F`NImpLen'.`price' `instrument' if `sample'
		gen esample = e(sample)
	}
	else {
		gen esample = 0
		replace esample = 1 if m`price'_price_ref_h`NImpLen'_r_SA < . & `sample'
	}
}

* Obtain sample length
sum dm if esample == 1
local firstMonth : di %tm = r(min)
local lastMonth : di %tm = r(max)
				
di "Sample: `firstMonth' - `lastMonth'"

preserve

foreach h of numlist 0/`NImpLen' {
	* Local projection step
	if "`level_or_cumulative'" == "level" {
		** Generate summed dependent variable (over horizon) ..
		* Sum up leads
		local prevH = `h' - 1
		if `h' == 0 gen `price'_h`h' = `price'
		* else gen `price'_h`h' = (F`h'.`price' + `price'_h`prevH'*`h')/(`h'+1)  //if inflation
		else gen `price'_h`h' = F`h'.`price' + `price'_h`prevH'

		local inflationToUse = "`price'_h`h'"
	}
	else { //if cumulative
		local inflationToUse = "m`price'_price_ref_h`h'_r_SA"
	}
					
	* Instrument is contemporaneous (this is a hack to ensure we can keep the same code for the regressions)(all it does is say "use the instrument as is")
	local instrumentToUse = "`instrument'"
	 
	newey `inflationToUse' `instrumentToUse' `controls' if esample == 1, lag(24)
	
	* Store results in matrix
	local rowNumber = `h' + 1
					
	mat b = e(b)
	mat V = e(V)
	mat se = sqrt(V[1,1])

	mat results[`rowNumber', 1] = `h'
	mat results[`rowNumber', 2] = b[1,1]
	mat results[`rowNumber', 3] = se[1,1]
}

* Store results in dataset
clear
svmat results, names(col)
gen ci_90L = b - 1.645 * se
gen ci_90H = b + 1.645 * se
gen ci_67L = b - 0.98 * se
gen ci_67H = b + 0.98 * se
gen ci_95L = b - 1.96 * se
gen ci_95H = b + 1.96 * se
* gen instrument = "`instrument'"
* gen price = "`price'"
				
* Plot
local conf "95"  //or "90" or "67"
twoway line ci_`conf'L b ci_`conf'H horizon, lcolor(black%30 maroon black%30) yline(0, lcolor(black%75) lpattern(.-)) legend(off) ytitle("%", size(large)) xtitle("months", size(large)) xlabel(, labsize(large)) ylabel(, labsize(large)) graphregion(color(white)) //note("`price_label' - `instrument_label' - Sample: `firstMonth' - `lastMonth' - 90% bands" "Controls: `controls_label'") 
graph export "$path_graphs/`price'_`instrument'_`controls_name'_`sample_str'_`conf'_newey.pdf", replace as(pdf)

restore

* log close
